Rarity models (with IUCN evaluated species)
Simple rarity model
if (!file.exists("./vet_rare_naturalized_10k.rds")) {
vet_rare_nat <- brm(rare | trials(vetted) ~ Naturalised_scaled + diversification_scaled + new_age_scaled + (1|family) + (1|family_name),
dat=data_vetted, family=binomial(),
iter=10000, thin=4, cov_ranef = list(family = phy_cov_vetted),
control=list(adapt_delta=0.95,max_treedepth=12), cores=4)
saveRDS(vet_rare_nat, "./vet_rare_naturalized_10k.rds")
} else { vet_rare_nat <- readRDS("./vet_rare_naturalized_10k.rds")}
outcome <- summary(vet_rare_nat$fit)
to_print <- outcome$summary[grep("^[sd,b]_*",rownames(outcome$summary)),]
to_print <- to_print[-1,]
names1 <- c("Number of naturalized species","Diversification rate","Family age","SD brownian effects","SD family specific effects")
vet_rare_nat_plot <- stanplot(vet_rare_nat, pars=rownames(to_print), type="intervals",
prob=0.8, prob_outer=0.95, point_est="mean") + scale_y_discrete(labels=rev(names1), limits=rev(rownames(to_print)))
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
vet_rare_nat_plot

# pdf("./plots_tables/vet_rare_nat_plot.pdf", width=5, height=4)
# vet_rare_nat_plot + ggtitle("")
# dev.off()
rownames(to_print) <- names1
kable(to_print, digits=2)
| Number of naturalized species |
-0.98 |
0.00 |
0.20 |
-1.38 |
-1.12 |
-0.98 |
-0.84 |
-0.57 |
4310.33 |
1.00 |
| Diversification rate |
1.67 |
0.01 |
0.35 |
0.99 |
1.44 |
1.66 |
1.90 |
2.34 |
4418.13 |
1.00 |
| Family age |
1.32 |
0.01 |
0.33 |
0.66 |
1.11 |
1.33 |
1.55 |
1.96 |
4159.95 |
1.00 |
| SD brownian effects |
0.98 |
0.01 |
0.20 |
0.60 |
0.84 |
0.97 |
1.12 |
1.37 |
1048.18 |
1.00 |
| SD family specific effects |
0.59 |
0.00 |
0.14 |
0.25 |
0.51 |
0.61 |
0.69 |
0.83 |
922.67 |
1.01 |
to_print <- to_print[,c("mean","sd","2.5%","50%","97.5%","n_eff","Rhat")]
print(xtable(to_print, type = "latex"), file = "./plots_tables/vet_rare_nat.tex")
##Posterior predictive checks
pp_vet_rare_nat <- brms::posterior_predict(vet_rare_nat)
log1 <- scale_x_continuous(trans="log1p")
ppc_dens_overlay(y=vet_rare_nat$data$rare, pp_vet_rare_nat[1:50, ]) + scale_x_continuous(trans="log1p")

# NRMSE
nrmse_vet_rare_nat <- apply(pp_vet_rare_nat, 1, function(x) nrmse(pred=x,obs=vet_rare_nat$data$rare, type="iq"))
mean(nrmse_vet_rare_nat)
## [1] 0.1650106
sd(nrmse_vet_rare_nat)
## [1] 0.01928937
# family-level variance due to brownian motion component
hyp <- paste("sd_family__Intercept^2 /", "(sd_family__Intercept^2 + sd_family_name__Intercept^2) = 0")
(hyp <- hypothesis(vet_rare_nat, hyp, class = NULL))
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (sd_family__Inter... = 0 0.71 0.16 0.36 0.96 NA
## Post.Prob Star
## 1 NA *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
Full rarity model (IUCN vetted species)
if (!file.exists("./vet_rare_10k.rds")) {
vet_rare <- brm(rare | trials(vetted) ~ diversification_scaled + new_age_scaled + woodiness + herbaceous + climate.sum + species.hybrids + hermaphroditic + dioecious + monoecious + asexual + fleshy.fruit + animal + (1|family) + (1|family_name),
dat=data_vetted, family=binomial(),
iter=10000, thin=4, cov_ranef = list(family = phy_cov_vetted),
control=list(adapt_delta=0.95,max_treedepth=13))
saveRDS(vet_rare, "./vet_rare_10k.rds")
} else { vet_rare <- readRDS("./vet_rare_10k.rds")}
# prior_summary(vet_rare, all=FALSE)
outcome <- summary(vet_rare$fit)
to_print <- outcome$summary[grep("^[sd,b]_*",rownames(outcome$summary)),]
to_print <- to_print[-1,]
names1 <- c("Diversification rate","Famly age","Woodiness","Herbaceous","Climate sum","Number of species hybrids","Hermaphroditic","Dioecious","Monoecious","Asexual","Fleshy fruit","Animal pollinated","SD brownian effects","SD family specific effects")
vet_rare_plot <- stanplot(vet_rare, pars=rownames(to_print), type="intervals",
prob=0.8, prob_outer=0.95, point_est="mean") + scale_y_discrete(labels=rev(names1), limits=rev(rownames(to_print))) + ggtitle("Proportion rare species")
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
vet_rare_plot

# pdf("./plots_tables/vet_rare_plot.pdf", width=5, height=4)
# vet_rare_plot + ggtitle("")
# dev.off()
rownames(to_print) <- names1
kable(to_print, digits=2)
| Diversification rate |
0.51 |
0.01 |
0.32 |
-0.13 |
0.30 |
0.51 |
0.72 |
1.13 |
3964.93 |
1.00 |
| Famly age |
0.48 |
0.00 |
0.29 |
-0.09 |
0.28 |
0.48 |
0.67 |
1.04 |
3359.11 |
1.00 |
| Woodiness |
0.61 |
0.00 |
0.13 |
0.36 |
0.52 |
0.61 |
0.70 |
0.87 |
4220.74 |
1.00 |
| Herbaceous |
-1.10 |
0.00 |
0.21 |
-1.51 |
-1.24 |
-1.10 |
-0.95 |
-0.67 |
2896.15 |
1.00 |
| Climate sum |
-0.25 |
0.00 |
0.09 |
-0.43 |
-0.31 |
-0.25 |
-0.19 |
-0.07 |
2404.65 |
1.00 |
| Number of species hybrids |
-0.18 |
0.00 |
0.13 |
-0.44 |
-0.27 |
-0.18 |
-0.09 |
0.07 |
3256.32 |
1.00 |
| Hermaphroditic |
0.12 |
0.01 |
0.22 |
-0.34 |
-0.02 |
0.13 |
0.28 |
0.56 |
1356.28 |
1.00 |
| Dioecious |
-0.13 |
0.00 |
0.15 |
-0.42 |
-0.23 |
-0.13 |
-0.03 |
0.17 |
4247.49 |
1.00 |
| Monoecious |
-0.13 |
0.00 |
0.15 |
-0.44 |
-0.24 |
-0.14 |
-0.03 |
0.17 |
4252.41 |
1.00 |
| Asexual |
0.05 |
0.00 |
0.16 |
-0.25 |
-0.05 |
0.05 |
0.16 |
0.37 |
4154.96 |
1.00 |
| Fleshy fruit |
0.06 |
0.00 |
0.14 |
-0.23 |
-0.04 |
0.05 |
0.15 |
0.34 |
4643.26 |
1.00 |
| Animal pollinated |
0.55 |
0.00 |
0.24 |
0.06 |
0.38 |
0.55 |
0.72 |
1.01 |
3453.40 |
1.00 |
| SD brownian effects |
0.47 |
0.01 |
0.27 |
0.03 |
0.25 |
0.46 |
0.66 |
1.02 |
455.58 |
1.01 |
| SD family specific effects |
0.75 |
0.00 |
0.11 |
0.49 |
0.69 |
0.76 |
0.82 |
0.92 |
685.94 |
1.00 |
to_print <- to_print[,c("mean","sd","2.5%","50%","97.5%","n_eff","Rhat")]
print(xtable(to_print, type = "latex"), file = "./plots_tables/vet_rare.tex")
##Posterior predictive checks
# first make posterior predictions
pp_rare <- brms::posterior_predict(vet_rare)
log1 <- scale_x_continuous(trans="log1p")
ppc_dens_overlay(y=vet_rare$data$rare, pp_rare[1:50, ]) + scale_x_continuous(trans="log1p")

# NRMSE
nrmse_vet <- apply(pp_rare, 1, function(x) nrmse(pred=x,obs=vet_rare$data$rare, type="iq"))
mean(nrmse_vet)
## [1] 0.164503
sd(nrmse_vet)
## [1] 0.01905102
# family-level variance due to brownian motion component
hyp <- paste("sd_family__Intercept^2 /", "(sd_family__Intercept^2 + sd_family_name__Intercept^2) = 0")
(hyp <- hypothesis(vet_rare, hyp, class = NULL))
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (sd_family__Inter... = 0 0.3 0.24 0 0.8 NA
## Post.Prob Star
## 1 NA *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# plot(hyp)
Naturalized model
if (!file.exists("./fit_nat_10k.rds")) {
fit_nat <- brm(Naturalised | trials(species) ~ diversification_scaled + new_age_scaled + woodiness + herbaceous + climate.sum + species.hybrids + hermaphroditic + dioecious + monoecious + asexual + fleshy.fruit + animal + (1|family) + (1|family_name),
dat=data, family=binomial(),
iter=10000, thin=4, cov_ranef = list(family = phy_cov),
control=list(adapt_delta=0.8, max_treedepth=11))
saveRDS(fit_nat, "./fit_nat_10k.rds")
} else { fit_nat <- readRDS("./fit_nat_10k.rds")}
outcome <- summary(fit_nat$fit)
to_print <- outcome$summary[grep("^[sd,b]_*",rownames(outcome$summary)),]
to_print <- to_print[-1,]
names1 <- c("Diversification rate","Famly age", "Woodiness","Herbaceous","Climate sum","Number of species hybrids","Hermaphroditic","Dioecious","Monoecious","Asexual","Fleshy fruit","Animal pollinated","SD brownian effects","SD family specific effects")
fit_nat_plot <- stanplot(fit_nat, pars=rownames(to_print), type="intervals",
prob=0.8, prob_outer=0.95, point_est="mean") + scale_y_discrete(labels=rev(names1), limits=rev(rownames(to_print)))+ ggtitle("Proportion of naturalized species")
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
fit_nat_plot

# pdf("./plots_tables/fit_nat_plot.pdf", width=5, height=4)
# fit_nat_plot + ggtitle("")
# dev.off()
rownames(to_print) <- names1
kable(to_print, digits=2)
| Diversification rate |
-1.32 |
0.00 |
0.27 |
-1.86 |
-1.51 |
-1.33 |
-1.14 |
-0.78 |
4163.57 |
1 |
| Famly age |
-0.94 |
0.00 |
0.30 |
-1.53 |
-1.14 |
-0.95 |
-0.75 |
-0.35 |
4640.06 |
1 |
| Woodiness |
-0.05 |
0.00 |
0.14 |
-0.32 |
-0.14 |
-0.04 |
0.05 |
0.22 |
4827.86 |
1 |
| Herbaceous |
0.92 |
0.00 |
0.23 |
0.47 |
0.77 |
0.91 |
1.07 |
1.38 |
4976.75 |
1 |
| Climate sum |
0.27 |
0.00 |
0.09 |
0.09 |
0.21 |
0.27 |
0.34 |
0.46 |
4835.06 |
1 |
| Number of species hybrids |
0.40 |
0.00 |
0.12 |
0.16 |
0.32 |
0.40 |
0.48 |
0.64 |
4905.79 |
1 |
| Hermaphroditic |
0.05 |
0.00 |
0.24 |
-0.42 |
-0.11 |
0.05 |
0.22 |
0.53 |
4234.42 |
1 |
| Dioecious |
-0.02 |
0.00 |
0.16 |
-0.33 |
-0.13 |
-0.02 |
0.09 |
0.31 |
4870.42 |
1 |
| Monoecious |
0.01 |
0.00 |
0.16 |
-0.31 |
-0.10 |
0.01 |
0.12 |
0.32 |
4747.97 |
1 |
| Asexual |
0.50 |
0.00 |
0.17 |
0.17 |
0.39 |
0.50 |
0.62 |
0.84 |
4797.08 |
1 |
| Fleshy fruit |
-0.13 |
0.00 |
0.14 |
-0.40 |
-0.23 |
-0.13 |
-0.04 |
0.14 |
5196.49 |
1 |
| Animal pollinated |
-0.49 |
0.00 |
0.28 |
-1.04 |
-0.68 |
-0.48 |
-0.30 |
0.07 |
4930.34 |
1 |
| SD brownian effects |
1.53 |
0.01 |
0.18 |
1.14 |
1.41 |
1.54 |
1.65 |
1.83 |
790.83 |
1 |
| SD family specific effects |
0.36 |
0.01 |
0.21 |
0.02 |
0.19 |
0.36 |
0.51 |
0.75 |
525.43 |
1 |
to_print <- to_print[,c("mean","sd","2.5%","50%","97.5%","n_eff","Rhat")]
print(xtable(to_print, type = "latex"), file = "./plots_tables/fit_nat.tex")
##Posterior predictive checks
pp_nat <- brms::posterior_predict(fit_nat)
log1 <- scale_x_continuous(trans="log1p")
ppc_dens_overlay(y=fit_nat$data$Naturalised, pp_nat[1:50, ]) + scale_x_continuous(trans="log1p")

# NRMSE
nrmse_nat <- apply(pp_nat, 1, function(x) nrmse(pred=x,obs=fit_nat$data$Naturalised, type="iq"))
mean(nrmse_nat)
## [1] 0.5494844
sd(nrmse_nat)
## [1] 0.0757893
# family-level variance due to brownian motion component
hyp <- paste("sd_family__Intercept^2 /", "(sd_family__Intercept^2 + sd_family_name__Intercept^2) = 0")
(hyp <- hypothesis(fit_nat, hyp, class = NULL))
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (sd_family__Inter... = 0 0.92 0.08 0.71 1 NA
## Post.Prob Star
## 1 NA *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# simplified threat model
outcome <- summary(vet_rare_nat$fit, prob=c(0.1,0.9))
to_print <- outcome$summary[grep("^[sd,b]_*",rownames(outcome$summary)),]
to_print <- to_print[-1,]
to_print <- data.frame(to_print)
names1 <- c("Number of naturalized species","Diversification rate","Family age","SD Brownian effects","SD family specific effects")
to_print$variable <- as.factor(names1)
p <- ggplot(to_print, aes(x=variable, y=mean)) + scale_x_discrete(limits = rev( levels(to_print$variable))) +
geom_pointrange(aes(ymin=X10., ymax=X90.)) + coord_flip() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
geom_hline(yintercept = 0) + xlab("") + ylab("") +
theme(legend.title = element_blank())
p

# ggsave("./plots_tables/threat_nat_forest.pdf", p , height=3,width=6)
# ggsave("./plots_tables/threat_nat_forest.png", p , height=3,width=6)
# combined Threat and Naturalized
threat_out <- summary(vet_rare$fit, prob=c(0.1,0.9))
threat_out <- threat_out$summary[grep("^[sd,b]_*",rownames(threat_out$summary)),]
threat_out <- threat_out[-1,]
threat_out <- data.frame(threat_out)
threat_out$group <- rev(letters[1:nrow(threat_out)])
threat_out$model <- "Threatened"
nat_out <- summary(fit_nat$fit, prob=c(0.1,0.9))
nat_out <- nat_out$summary[grep("^[sd,b]_*",rownames(nat_out$summary)),]
nat_out <- nat_out[-1,]
nat_out <- data.frame(nat_out)
nat_out$group <- rev(letters[1:nrow(nat_out)])
nat_out$model <- "Naturalized"
labels=c("Diversification rate","Family age","Variation in growth form","Herbaceous","Climate range variability","Species hybrids","Hermaphroditic","Dioecious","Monoecious","Asexual seed production","Fleshy fruits","Animal pollination","SD Brownian effects","SD family specific effects")
threat_out$variable <- labels
nat_out$variable <- labels
coefs <- rbind(threat_out, nat_out)
coefs$overlap <- c( "n","n","n","n","n","n","y","y","y","y","y","n","n","n",
"n","n","y","n","n","n","y","y","y","n","y","n","n","n")
coefs$opposite <- c("y","y","n","y","y","y","n","n","n","n","n","y","n","n",
"y","y","n","y","y","y","n","n","n","n","n","y","n","n")
coefs$both <- c("ny","ny","nn","ny","ny","ny","yn","yn","yn","yn","yn","ny","nn","nn",
"ny","ny","yn","ny","ny","ny","yn","yn","yn","nn","yn","ny","nn","nn")
labels<-rev(labels)
cols<-c("#88b7c5","#a36666")
p2 <- ggplot(coefs, aes(x=group, y=mean, color=model)) + scale_shape_manual(values=c(17,15,17)) + scale_alpha_manual(values=c(1,1,0.4)) +
geom_pointrange(aes(ymin=X10., ymax=X90.,shape=both,alpha=both)) + coord_flip() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black")) +
geom_hline(yintercept = 0) + xlab("") + ylab("") + theme(legend.position = "none") +
theme(legend.title = element_blank()) + scale_x_discrete(labels=labels) + guides(shape = FALSE, size = FALSE) + scale_colour_manual(values=cols)
p2

# ggsave("./plots_tables/combined_forest.pdf", p2 , height=5,width=6)
# ggsave("./plots_tables/combined_forest.png", p2 , height=5,width=6)


Sensitivity analyses
Full rarity model: no family level effects
if (!file.exists("./vet_rare_nofam_10k.rds")) {
vet_rare_nofam <- brm(rare | trials(vetted) ~ diversification_scaled + new_age_scaled + woodiness + herbaceous + climate.sum + species.hybrids + hermaphroditic + dioecious + monoecious + asexual + fleshy.fruit + animal,
dat=data_vetted, family=binomial(),
iter=10000, thin=4,
control=list(adapt_delta=0.95,max_treedepth=13))
saveRDS(vet_rare_nofam, "./vet_rare_nofam_10k.rds")
} else { vet_rare_nofam <- readRDS("./vet_rare_nofam_10k.rds")}
outcome <- summary(vet_rare_nofam$fit)
to_print <- outcome$summary[grep("^[sd,b]_*",rownames(outcome$summary)),]
to_print <- to_print[-1,]
names1 <- c("Diversification rate","Famly age", "Woodiness","Herbaceous","Climate sum","Number of species hybrids","Hermaphroditic","Dioecious","Monoecious","Asexual","Fleshy fruit","Animal pollinated")
vet_rare_nofam_plot <- stanplot(vet_rare_nofam, pars=rownames(to_print), type="intervals",
prob=0.8, prob_outer=0.95, point_est="mean") + scale_y_discrete(labels=rev(names1), limits=rev(rownames(to_print))) + ggtitle("Proportion rare species")
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
vet_rare_nofam_plot

# pdf("./plots_tables/vet_rare_nofam_plot.pdf", width=5, height=4)
# vet_rare_nofam_plot + ggtitle("")
# dev.off()
rownames(to_print) <- names1
kable(to_print, digits=2)
| Diversification rate |
0.22 |
0 |
0.09 |
0.05 |
0.16 |
0.22 |
0.28 |
0.40 |
4714.34 |
1 |
| Famly age |
0.24 |
0 |
0.08 |
0.08 |
0.18 |
0.24 |
0.29 |
0.40 |
4714.28 |
1 |
| Woodiness |
0.48 |
0 |
0.03 |
0.43 |
0.46 |
0.48 |
0.50 |
0.54 |
4821.16 |
1 |
| Herbaceous |
-0.73 |
0 |
0.05 |
-0.83 |
-0.76 |
-0.73 |
-0.69 |
-0.63 |
5241.34 |
1 |
| Climate sum |
-0.35 |
0 |
0.02 |
-0.39 |
-0.36 |
-0.35 |
-0.33 |
-0.30 |
4903.43 |
1 |
| Number of species hybrids |
-0.15 |
0 |
0.04 |
-0.24 |
-0.18 |
-0.15 |
-0.12 |
-0.07 |
4562.55 |
1 |
| Hermaphroditic |
0.18 |
0 |
0.06 |
0.07 |
0.14 |
0.18 |
0.22 |
0.30 |
4969.73 |
1 |
| Dioecious |
0.14 |
0 |
0.03 |
0.07 |
0.11 |
0.14 |
0.16 |
0.20 |
4822.02 |
1 |
| Monoecious |
-0.25 |
0 |
0.04 |
-0.32 |
-0.27 |
-0.25 |
-0.22 |
-0.17 |
4751.92 |
1 |
| Asexual |
0.17 |
0 |
0.03 |
0.11 |
0.15 |
0.17 |
0.19 |
0.24 |
5007.34 |
1 |
| Fleshy fruit |
-0.24 |
0 |
0.02 |
-0.28 |
-0.25 |
-0.23 |
-0.22 |
-0.20 |
5019.82 |
1 |
| Animal pollinated |
0.93 |
0 |
0.06 |
0.81 |
0.88 |
0.92 |
0.97 |
1.04 |
4911.03 |
1 |
to_print <- to_print[,c("mean","sd","2.5%","50%","97.5%","n_eff","Rhat")]
print(xtable(to_print, type = "latex"), file = "./plots_tables/vet_rare_nofam.tex")
##Posterior predictive checks
# first make posterior predictions
pp_rare_nofam <- brms::posterior_predict(vet_rare_nofam)
log1 <- scale_x_continuous(trans="log1p")
ppc_dens_overlay(y=vet_rare_nofam$data$rare, pp_rare_nofam[1:50, ]) + scale_x_continuous(trans="log1p")

# NRMSE
nrmse_rare_nofam <- apply(pp_rare_nofam, 1, function(x) nrmse(pred=x,obs=vet_rare_nofam$data$rare, type="iq"))
mean(nrmse_rare_nofam)
## [1] 0.669739
sd(nrmse_rare_nofam)
## [1] 0.03574458
Full rarity model: only non-brownian family effects
if (!file.exists("./vet_rare_nophylo_10k.rds")) {
vet_rare_nophylo <- brm(rare | trials(vetted) ~ diversification_scaled + new_age_scaled + woodiness + herbaceous + climate.sum + species.hybrids + hermaphroditic + dioecious + monoecious + asexual + fleshy.fruit + animal + (1|family_name),
dat=data_vetted, family=binomial(),
iter=10000, thin=4,
control=list(adapt_delta=0.95,max_treedepth=13))
saveRDS(vet_rare_nophylo, "./vet_rare_nophylo_10k.rds")
} else { vet_rare_nophylo <- readRDS("./vet_rare_nophylo_10k.rds")}
outcome <- summary(vet_rare_nophylo$fit)
to_print <- outcome$summary[grep("^[sd,b]_*",rownames(outcome$summary)),]
to_print <- to_print[-1,]
names1 <- c("Diversification rate","Famly age", "Woodiness","Herbaceous","Climate sum","Number of species hybrids","Hermaphroditic","Dioecious","Monoecious","Asexual","Fleshy fruit","Animal pollinated","SD family specific effects")
vet_rare_nophylo_plot <- stanplot(vet_rare_nophylo, pars=rownames(to_print), type="intervals",
prob=0.8, prob_outer=0.95, point_est="mean") + scale_y_discrete(labels=rev(names1), limits=rev(rownames(to_print))) + ggtitle("Proportion rare species")
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
vet_rare_nophylo_plot

# pdf("./plots_tables/vet_rare_nophylo_plot.pdf", width=5, height=4)
# vet_rare_nophylo_plot + ggtitle("")
# dev.off()
rownames(to_print) <- names1
kable(to_print, digits=2)
| Diversification rate |
0.50 |
0.01 |
0.31 |
-0.10 |
0.29 |
0.50 |
0.71 |
1.11 |
3082.51 |
1 |
| Famly age |
0.55 |
0.00 |
0.27 |
0.01 |
0.36 |
0.55 |
0.72 |
1.08 |
3657.42 |
1 |
| Woodiness |
0.63 |
0.00 |
0.12 |
0.39 |
0.55 |
0.63 |
0.71 |
0.86 |
2945.50 |
1 |
| Herbaceous |
-1.18 |
0.00 |
0.18 |
-1.55 |
-1.30 |
-1.18 |
-1.06 |
-0.83 |
3948.81 |
1 |
| Climate sum |
-0.26 |
0.00 |
0.09 |
-0.44 |
-0.32 |
-0.26 |
-0.20 |
-0.08 |
3562.52 |
1 |
| Number of species hybrids |
-0.18 |
0.00 |
0.13 |
-0.44 |
-0.26 |
-0.18 |
-0.09 |
0.07 |
4150.29 |
1 |
| Hermaphroditic |
0.14 |
0.00 |
0.21 |
-0.27 |
-0.01 |
0.14 |
0.28 |
0.57 |
4096.12 |
1 |
| Dioecious |
-0.13 |
0.00 |
0.15 |
-0.42 |
-0.24 |
-0.14 |
-0.03 |
0.17 |
3526.58 |
1 |
| Monoecious |
-0.13 |
0.00 |
0.15 |
-0.43 |
-0.23 |
-0.13 |
-0.03 |
0.17 |
3403.72 |
1 |
| Asexual |
0.07 |
0.00 |
0.15 |
-0.24 |
-0.03 |
0.07 |
0.18 |
0.37 |
2911.12 |
1 |
| Fleshy fruit |
0.07 |
0.00 |
0.15 |
-0.22 |
-0.03 |
0.07 |
0.17 |
0.36 |
4293.12 |
1 |
| Animal pollinated |
0.60 |
0.00 |
0.23 |
0.16 |
0.44 |
0.60 |
0.76 |
1.07 |
4069.20 |
1 |
| SD family specific effects |
0.83 |
0.00 |
0.06 |
0.72 |
0.79 |
0.83 |
0.87 |
0.96 |
3781.04 |
1 |
to_print <- to_print[,c("mean","sd","2.5%","50%","97.5%","n_eff","Rhat")]
print(xtable(to_print, type = "latex"), file = "./plots_tables/vet_rare_nophylo.tex")
##Posterior predictive checks
pp_rare_nophylo <- brms::posterior_predict(vet_rare_nophylo)
log1 <- scale_x_continuous(trans="log1p")
ppc_dens_overlay(y=vet_rare_nophylo$data$rare, pp_rare_nophylo[1:50, ]) + scale_x_continuous(trans="log1p")

# NRMSE
nrmse_rare_nophylo <- apply(pp_rare_nophylo, 1, function(x) nrmse(pred=x,obs=vet_rare_nophylo$data$rare, type="iq"))
mean(nrmse_rare_nophylo)
## [1] 0.1648308
sd(nrmse_rare_nophylo)
## [1] 0.01937899
Full rarity model: only brownian family effects
if (!file.exists("./vet_rare_phylo_10k.rds")) {
vet_rare_phylo <- brm(rare | trials(vetted) ~ diversification_scaled + new_age_scaled + woodiness + herbaceous + climate.sum + species.hybrids + hermaphroditic + dioecious + monoecious + asexual + fleshy.fruit + animal+ (1|family),
dat=data_vetted, family=binomial(),
iter=10000, thin=4, cov_ranef = list(family = phy_cov),
control=list(adapt_delta=0.95,max_treedepth=13))
saveRDS(vet_rare_phylo, "./vet_rare_phylo_10k.rds")
} else { vet_rare_phylo <- readRDS("./vet_rare_phylo_10k.rds")}
outcome <- summary(vet_rare_phylo$fit)
to_print <- outcome$summary[grep("^[sd,b]_*",rownames(outcome$summary)),]
to_print <- to_print[-1,]
names1 <- c("Diversification rate","Famly age", "Woodiness","Herbaceous","Climate sum","Number of species hybrids","Hermaphroditic","Dioecious","Monoecious","Asexual","Fleshy fruit","Animal pollinated","SD brownian effects")
vet_rare_phylo_plot <- stanplot(vet_rare_phylo, pars=rownames(to_print), type="intervals",
prob=0.8, prob_outer=0.95, point_est="mean") + scale_y_discrete(labels=rev(names1), limits=rev(rownames(to_print))) + ggtitle("Proportion rare species")
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
vet_rare_phylo_plot

# pdf("./plots_tables/vet_rare_phylo_plot.pdf", width=5, height=4)
# vet_rare_phylo_plot + ggtitle("")
# dev.off()
rownames(to_print) <- names1
kable(to_print, digits=2)
| Diversification rate |
0.61 |
0.01 |
0.32 |
-0.01 |
0.39 |
0.62 |
0.83 |
1.23 |
3839.79 |
1 |
| Famly age |
0.42 |
0.01 |
0.33 |
-0.23 |
0.20 |
0.42 |
0.64 |
1.04 |
3398.65 |
1 |
| Woodiness |
0.57 |
0.00 |
0.13 |
0.31 |
0.48 |
0.57 |
0.66 |
0.84 |
3084.68 |
1 |
| Herbaceous |
-0.96 |
0.00 |
0.22 |
-1.40 |
-1.11 |
-0.96 |
-0.82 |
-0.53 |
2776.99 |
1 |
| Climate sum |
-0.22 |
0.00 |
0.09 |
-0.40 |
-0.28 |
-0.22 |
-0.16 |
-0.05 |
3085.21 |
1 |
| Number of species hybrids |
-0.19 |
0.00 |
0.12 |
-0.42 |
-0.27 |
-0.19 |
-0.10 |
0.05 |
4445.11 |
1 |
| Hermaphroditic |
0.06 |
0.00 |
0.22 |
-0.35 |
-0.09 |
0.06 |
0.22 |
0.49 |
3500.32 |
1 |
| Dioecious |
-0.16 |
0.00 |
0.15 |
-0.45 |
-0.26 |
-0.16 |
-0.06 |
0.13 |
3156.57 |
1 |
| Monoecious |
-0.09 |
0.00 |
0.15 |
-0.38 |
-0.19 |
-0.09 |
0.01 |
0.20 |
3288.50 |
1 |
| Asexual |
0.02 |
0.00 |
0.16 |
-0.30 |
-0.09 |
0.02 |
0.13 |
0.34 |
3022.65 |
1 |
| Fleshy fruit |
-0.07 |
0.00 |
0.12 |
-0.30 |
-0.15 |
-0.07 |
0.01 |
0.17 |
4274.28 |
1 |
| Animal pollinated |
0.44 |
0.00 |
0.25 |
-0.04 |
0.27 |
0.45 |
0.62 |
0.92 |
3181.12 |
1 |
| SD brownian effects |
1.31 |
0.00 |
0.10 |
1.13 |
1.24 |
1.31 |
1.37 |
1.51 |
3249.90 |
1 |
to_print <- to_print[,c("mean","sd","2.5%","50%","97.5%","n_eff","Rhat")]
print(xtable(to_print, type = "latex"), file = "./plots_tables/vet_rare_phylo.tex")
##Posterior predictive checks
pp_rare_phylo <- brms::posterior_predict(vet_rare_phylo)
log1 <- scale_x_continuous(trans="log1p")
ppc_dens_overlay(y=vet_rare_phylo$data$rare, pp_rare_phylo[1:50, ]) + scale_x_continuous(trans="log1p")

# NRMSE
nrmse_rare_phylo <- apply(pp_rare_phylo, 1, function(x) nrmse(pred=x,obs=vet_rare_phylo$data$rare, type="iq"))
mean(nrmse_rare_phylo)
## [1] 0.1651364
sd(nrmse_rare_phylo)
## [1] 0.01911785
Full rarity model (total species richness - not vetted by IUCN)
if (!file.exists("./nonvet_rare_10k.rds")) {
nonvet_rare <- brm(rare | trials(species) ~ diversification_scaled + new_age_scaled + woodiness + herbaceous + climate.sum + species.hybrids + hermaphroditic + dioecious + monoecious + asexual + fleshy.fruit + animal + (1|family) + (1|family_name),
dat=data, family=binomial(),
iter=10000, thin=4, cov_ranef = list(family = phy_cov),
control=list(adapt_delta=0.95,max_treedepth=13))
saveRDS(nonvet_rare, "./nonvet_rare_10k.rds")
} else { nonvet_rare <- readRDS("./nonvet_rare_10k.rds")}
outcome <- summary(nonvet_rare$fit)
to_print <- outcome$summary[grep("^[sd,b]_*",rownames(outcome$summary)),]
to_print <- to_print[-1,]
names1 <- c("Diversification rate","Famly age", "Woodiness","Herbaceous","Climate sum","Number of species hybrids","Hermaphroditic","Dioecious","Monoecious","Asexual","Fleshy fruit","Animal pollinated","SD brownian effects","SD family specific effects")
nonvet_rare_plot <- stanplot(nonvet_rare, pars=rownames(to_print), type="intervals",
prob=0.8, prob_outer=0.95, point_est="mean") + scale_y_discrete(labels=rev(names1), limits=rev(rownames(to_print))) + ggtitle("Proportion rare species")
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
nonvet_rare_plot

# pdf("./plots_tables/nonvet_rare_plot.pdf", width=5, height=4)
# nonvet_rare_plot + ggtitle("")
# dev.off()
rownames(to_print) <- names1
kable(to_print, digits=2)
| Diversification rate |
0.84 |
0.01 |
0.36 |
0.15 |
0.59 |
0.84 |
1.08 |
1.53 |
3754.44 |
1.00 |
| Famly age |
0.76 |
0.00 |
0.33 |
0.14 |
0.54 |
0.76 |
0.98 |
1.41 |
4772.19 |
1.00 |
| Woodiness |
-0.38 |
0.00 |
0.16 |
-0.70 |
-0.49 |
-0.38 |
-0.27 |
-0.06 |
3778.61 |
1.00 |
| Herbaceous |
-0.65 |
0.00 |
0.26 |
-1.16 |
-0.83 |
-0.65 |
-0.48 |
-0.15 |
3478.66 |
1.00 |
| Climate sum |
-0.21 |
0.00 |
0.11 |
-0.43 |
-0.29 |
-0.21 |
-0.14 |
0.00 |
4541.01 |
1.00 |
| Number of species hybrids |
0.04 |
0.00 |
0.17 |
-0.28 |
-0.07 |
0.04 |
0.15 |
0.36 |
4212.00 |
1.00 |
| Hermaphroditic |
-0.18 |
0.00 |
0.26 |
-0.69 |
-0.36 |
-0.17 |
0.00 |
0.34 |
4825.84 |
1.00 |
| Dioecious |
-0.24 |
0.00 |
0.19 |
-0.61 |
-0.37 |
-0.24 |
-0.12 |
0.13 |
4515.76 |
1.00 |
| Monoecious |
0.07 |
0.00 |
0.18 |
-0.29 |
-0.05 |
0.07 |
0.20 |
0.43 |
4732.96 |
1.00 |
| Asexual |
0.12 |
0.00 |
0.21 |
-0.30 |
-0.02 |
0.12 |
0.26 |
0.52 |
4056.46 |
1.00 |
| Fleshy fruit |
0.31 |
0.00 |
0.18 |
-0.03 |
0.19 |
0.32 |
0.43 |
0.66 |
4623.42 |
1.00 |
| Animal pollinated |
0.18 |
0.01 |
0.33 |
-0.46 |
-0.04 |
0.18 |
0.39 |
0.83 |
3478.39 |
1.00 |
| SD brownian effects |
0.96 |
0.01 |
0.33 |
0.25 |
0.74 |
0.97 |
1.18 |
1.61 |
659.00 |
1.01 |
| SD family specific effects |
1.03 |
0.01 |
0.16 |
0.69 |
0.94 |
1.05 |
1.14 |
1.30 |
924.73 |
1.00 |
to_print <- to_print[,c("mean","sd","2.5%","50%","97.5%","n_eff","Rhat")]
print(xtable(to_print, type = "latex"), file = "./plots_tables/nonvet_rare.tex")
##Posterior predictive checks
pp_rare <- brms::posterior_predict(nonvet_rare)
log1 <- scale_x_continuous(trans="log1p")
ppc_dens_overlay(y=nonvet_rare$data$rare, pp_rare[1:50, ]) + scale_x_continuous(trans="log1p")

# NRMSE
nrmse_nonvet <- apply(pp_rare, 1, function(x) nrmse(pred=x,obs=nonvet_rare$data$rare, type="iq"))
mean(nrmse_nonvet)
## [1] 0.4346414
sd(nrmse_nonvet)
## [1] 0.04654522
# family-level variance due to brownian motion component
hyp <- paste("sd_family__Intercept^2 /", "(sd_family__Intercept^2 + sd_family_name__Intercept^2) = 0")
(hyp <- hypothesis(nonvet_rare, hyp, class = NULL))
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (sd_family__Inter... = 0 0.45 0.21 0.04 0.84 NA
## Post.Prob Star
## 1 NA *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
Invasive model
if (!file.exists("./fit_inv_10k.rds")) {
fit_inv <- brm(Invasive | trials(species) ~ diversification_scaled + new_age_scaled + woodiness + herbaceous + climate.sum + species.hybrids + hermaphroditic + dioecious + monoecious + asexual + fleshy.fruit + animal + (1|family) + (1|family_name),
dat=data, family=binomial(),
iter=10000, thin=4, cov_ranef = list(family = phy_cov),
control=list(adapt_delta=0.8, max_treedepth=11))
saveRDS(fit_inv, "./fit_inv_10k.rds")
} else { fit_inv <- readRDS("./fit_inv_10k.rds")}
outcome <- summary(fit_inv$fit)
to_print <- outcome$summary[grep("^[sd,b]_*",rownames(outcome$summary)),]
to_print <- to_print[-1,]
names1 <- c("Diversification rate","Famly age", "Woodiness","Herbaceous","Climate sum","Number of species hybrids","Hermaphroditic","Dioecious","Monoecious","Asexual","Fleshy fruit","Animal pollinated","SD brownian effects","SD family specific effects")
fit_inv_plot <- stanplot(fit_inv, pars=rownames(to_print), type="intervals",
prob=0.8, prob_outer=0.95, point_est="mean") + scale_y_discrete(labels=rev(names1), limits=rev(rownames(to_print)))+ ggtitle("Proportion of invasive species")
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
fit_inv_plot

# pdf("./plots_tables/fit_inv_plot.pdf", width=5, height=4)
# fit_inv_plot + ggtitle("")
# dev.off()
rownames(to_print) <- names1
kable(to_print, digits=2)
| Diversification rate |
-0.92 |
0.01 |
0.41 |
-1.69 |
-1.21 |
-0.94 |
-0.65 |
-0.11 |
4863.28 |
1.00 |
| Famly age |
-1.19 |
0.01 |
0.41 |
-1.99 |
-1.46 |
-1.19 |
-0.91 |
-0.39 |
4350.25 |
1.00 |
| Woodiness |
0.24 |
0.00 |
0.18 |
-0.11 |
0.12 |
0.24 |
0.35 |
0.59 |
4942.94 |
1.00 |
| Herbaceous |
0.33 |
0.00 |
0.31 |
-0.28 |
0.12 |
0.33 |
0.54 |
0.93 |
4609.73 |
1.00 |
| Climate sum |
0.05 |
0.00 |
0.13 |
-0.20 |
-0.04 |
0.05 |
0.13 |
0.31 |
4200.42 |
1.00 |
| Number of species hybrids |
0.49 |
0.00 |
0.16 |
0.16 |
0.38 |
0.49 |
0.59 |
0.80 |
4664.47 |
1.00 |
| Hermaphroditic |
-0.09 |
0.00 |
0.32 |
-0.71 |
-0.30 |
-0.09 |
0.13 |
0.54 |
5165.25 |
1.00 |
| Dioecious |
0.01 |
0.00 |
0.20 |
-0.40 |
-0.13 |
0.01 |
0.15 |
0.40 |
4781.75 |
1.00 |
| Monoecious |
0.06 |
0.00 |
0.21 |
-0.36 |
-0.08 |
0.06 |
0.19 |
0.46 |
4191.75 |
1.00 |
| Asexual |
0.42 |
0.00 |
0.21 |
0.00 |
0.28 |
0.42 |
0.56 |
0.83 |
4433.06 |
1.00 |
| Fleshy fruit |
0.00 |
0.00 |
0.18 |
-0.36 |
-0.12 |
0.00 |
0.13 |
0.36 |
4760.63 |
1.00 |
| Animal pollinated |
-0.22 |
0.01 |
0.36 |
-0.91 |
-0.46 |
-0.21 |
0.02 |
0.49 |
4701.72 |
1.00 |
| SD brownian effects |
1.09 |
0.02 |
0.44 |
0.12 |
0.81 |
1.15 |
1.43 |
1.77 |
676.64 |
1.01 |
| SD family specific effects |
0.66 |
0.01 |
0.29 |
0.05 |
0.47 |
0.70 |
0.88 |
1.12 |
695.02 |
1.01 |
to_print <- to_print[,c("mean","sd","2.5%","50%","97.5%","n_eff","Rhat")]
print(xtable(to_print, type = "latex"), file = "./plots_tables/invasive.tex")
##Posterior predictive checks
pp_inv <- brms::posterior_predict(fit_inv)
log1 <- scale_x_continuous(trans="log1p")
ppc_dens_overlay(y=fit_inv$data$Invasive, pp_inv[1:50, ]) + scale_x_continuous(trans="log1p")

# NRMSE
nrmse_inv <- apply(pp_inv, 1, function(x) nrmse(pred=x,obs=fit_inv$data$Invasive, type="iq"))
mean(nrmse_inv)
## [1] 1.445904
sd(nrmse_inv)
## [1] 0.2081148
# family-level variance due to brownian motion component
hyp <- paste("sd_family__Intercept^2 /", "(sd_family__Intercept^2 + sd_family_name__Intercept^2) = 0")
(hyp <- hypothesis(fit_inv, hyp, class = NULL))
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (sd_family__Inter... = 0 0.65 0.3 0.01 1 NA
## Post.Prob Star
## 1 NA *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.